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A new averaging method linking discrete to continuum variables of granular materials is developed 
and used to derive average balance equations. Its novelty lies in the choice of the decomposition 
between mean values and fluctuations of properties which takes into account the effect of gradients. 
Thanks to a local homogeneity hypothesis, whose validity is discussed, simplified balance equations 
are obtained. This original approach solves the problem of dependence of some variables on the size of 
the averaging domain obtained in previous approaches which can lead to huge relative errors (several 
hundred percentages). It also clearly separates affine and nonaffine fields in the balance equations. 
The resulting energy cascade picture is discussed, with a particular focus on unidirectional steady 
and fully developed flows for which it appears that the contact terms are dissipated locally unlike the 
kinetic terms which contribute to a nonlocal balance. Application of the method is demonstrated 
in the determination of the macroscopic properties such as volume fraction, velocity, stress, and 
energy of a simple shear flow, where the discrete results are generated by means of discrete particle 
simulation. 

PACS numbers: 45.70.-n, 47.57.Gc, 83.10.Ff 
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I. INTRODUCTION 


c3 The behavior of granular materials under deformation and 
Qow is often studied by means of numerical methods that 
is olve the equations of motion for every single body and ac¬ 
count for the interactions between particles through appro¬ 
priate models (see, e.g., Refs. [IH3|). These methods, com- 
i^ponly grouped under the name of discrete element methods 
L —(DEM) (j-jg, can be extremely useful for simulating discrete 
media with the purpose of establishing micro-macro relations, 
l^and in this perspective they need consistent averaging proce¬ 
ed ures [zHU- An averaging method computes continuum-like 
f^yariables (e.g., a velocity field) starting from their discrete, 
•^jparticle-scale counterparts (e.g., particle velocities). When 
OAsing such an approach, two issues have to be considered with 
^^itmost attention: the representativity of the average and the 
(^-physical meaning of the obtained estimates. We can say that 
(^^m average performed at a given point is representative if a 
lOufhciently large domain in space and time can be defined, 
^vhere particles share the same properties. Such a volume, 
^usually called “representative volume element,” may be de- 
UTjined with respect to only some variables. 

» i An important step in the development of averaging tech- 
Cdiiques was made by Babic [7]. In that work the author devel¬ 
oped a general framework for weighted space-time averages, 
applicable to a wide range of conditions, and giving a large 
set of self-consistent continuum balance equations. In Babic’s 
method, the definition of a representative volume is intrin¬ 
sically that of a zone where affine velocity fields are locally 
uniform, i.e., no gradients occur. The assumption that such 
a volume can be defined is called by Babic the “continuum 
assumption.” Now, the lack of scale separation typical of gran¬ 
ular flows (the scale of spatial variation of variables has the 
same order of magnitude of the particle size) implies that in 
a three-dimensional (3D) flow we cannot in practice define 
such a volume. As a consequence, the physical meaning of 
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averages obtained using Babic’s approach in a 3D flow is gen¬ 
erally questionable. At least some average estimates will then 
probably depend on the size of the averaging volume 10. 

We can see a simple example of this by considering a uni¬ 
form shear flow; if all particles follow the mean flow with no 
frustration, we will tend to say that fluctuations are negligi¬ 
ble. But this depends on the definition of fluctuation. If we 
calculate the fluctuational kinetic energy e T [which is related 
to granular temperature 0; see also Eq. (fl5l) ] by considering 
fluctuations as Babic does, we find e T ~ (7 D) 2 , where D is 
the averaging scale. Thus the fluctuations depend on the av¬ 
eraging scale which casts doubt on such a method. Also the 
kinetic contribution to the stress tensor will suffer the same 
problem. Thus, as we will show below, such an approach may 
lead to errors that can reach several hundred percentages. 

In this work we focus on this problem, with a particular 
attention on the mechanical energy equation. We discuss both 
the issue of dependence on the size of the averaging domain 
and the energy cascade picture associated with the averaging 
method. 

It is clear that the dependence of some average variables 
on the averaging domain size is intrinsic to Babic’s definition 
of fluctuations. In the following we show how a different fluc¬ 
tuation decomposition taking into account gradients yields a 
new derivation of continuum balance equations. These bal¬ 
ance equations have the merit of clearly separating affine and 
nonaffine fields and identifying the terms responsible for aver¬ 
aging scale dependence. With this decomposition it is possible 
to attribute a more sound physical meaning to the terms ap¬ 
pearing in the final balance equations. We show how these 
balance equations may be greatly simplified by a “local homo¬ 
geneity assumption,” which is strongly related to the fluctua¬ 
tion decomposition itself and which yields less strict require¬ 
ments than Babic’s “continuum assumption.” Applying this 
method to a simple class of flows, we find that two separate 
paths for fluctuating energy dissipation are implied: while ki¬ 
netic stress power enters a nonlocal balance, the contact stress 
power seems to be dissipated locally. This result, as well as 
the robustness of the method, is verified with discrete clement 
simulations of a simple shear flow. 



II. A NEW DERIVATION OF BALANCE 
EQUATIONS 


B. Continuum balances 
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A. Weighted average, fluctuation decomposition and 
local homogeneity assumption 

Given a particle property ip p (?tp,t'), Babic defined the av¬ 
eraged property at point x and time t as a weighted space- 
time average: 

/ OO 

y w(x p — x, t! — t)m p ip p dt ', (1) 

-oo p 

where w(x p — x, t' — t) = is a normalized weighting func¬ 
tion, p the average density, and m p the mass of a grain. The 
fluctuation was simply defined by Babic with respect to the 
center of the averaging volume: 


The local homogeneity hypothesis leaves the continuity 
equation unchanged: 

J^ + V-pv = 0, (5) 

where average density and velocity are defined by: 


/ OO 

w p m p dt ', 

-oo p 


pv = £ w p m p v p dt'. 

J-OC p 


( 6 ) 

( 7 ) 


^p = i) p - V’(x). (2) 

This may seem harmless, but in fact— as shown in 
Ref- 0—the dependence of some terms of the continuum av¬ 
eraged equations on the size of the averaging domain is strictly 
related to this choice for the fluctuation decomposition. 

We propose to define the fluctuation not with respect to 
the average value at point x (as was done in Ref. [7j), but 
with respect to the particle center, 

i>p=i> p -'${x p ). (3) 

In the following we will redevelop continuum-averaged equa¬ 
tions in the spirit of Babic’s approach, but using this new 
fluctuation decomposition. 

In order to simplify the result, we will make a constitutive 
assumption. We suppose that a scale exists where the 
gradients of ^ are smooth and that if this scale is the 
averaging scale intrinsic to the weighting function w pi we can 
approximate Vi/> as a constant near the averaging point. We 
call this assumption the local homogeneity assumption. To¬ 
gether with this assumption we also assume a homogeneous 
distribution of particle centers near the averaging point and 
the decorrelation of positions and velocities. These three 
assumptions are strictly related. 

While the absence of scale separation usually prevents iden¬ 
tification of a scale where gradients are zero, such a first-order 
approximation is reasonable and—as we will see—very fruit¬ 
ful. 

The local homogeneity assumption, involving gradients in¬ 
stead of variables, is less strict than Babic’s continuum hy¬ 
pothesis and therefore more likely to be valid under a proper 
choice of the weighting function. Based on this assumption, 
we can therefore approximate tH x p) by 

^(xp) = tA(x) + (x p -x)-V^(x). (4) 

Following Babic’s approach, it is possible to develop a 
generic balance equation for the particle property ip(’K p ). De¬ 
tails are given in Appendix A; here we show the main findings, 
with a particular regard to the continuity, momentum, and 
translational energy equations. When performing the deriva¬ 
tion, many new terms appear (with respect to Babic’s treat¬ 
ment), but most of them vanish due to the local homogeneity 
hypothesis. 


As for the momentum balance, we end up with the classical 
equation: 

d 

—pv + V • pvv = V • T + pg, (8) 

where the total stress tensor is composed by three contribu¬ 
tions: T = T c + T fc + T 7 . The first two terms are the same 
as in Babic. The first contribution is due to contact forces 
(both long-lasting contacts and collisions) and is given by 

/ OO 

££< Wp^', ( 9 ) 

'°° p q>p 

where w pq is a weight function related to the fraction of the 
branch vector joining the two particles p and q lying within the 
averaging function. Due to the new decomposition, Babic’s 
kinetic stress tensor is decomposed in a truly kinetic part, 

/ OO 

p m p v p v p dt', (10) 

-OO p 

and a new term, 

T 7 = —p(D • Vv)(D • Vv), (11) 

which contains all the dependence on averaging domain size 

through the vector D. It is useful to recall that the emergence 
of a part of the stress tensor depending on velocity gradients 
has nothing to do with constitutive relations but is the joint 
product of coarse graining and scale dependence. Due to their 
definition, both T 7 and T fc are symmetric. Evidently, in ab¬ 
sence of velocity gradients, T 7 will vanish; in that case, the 
kinetic stress tensor as defined by Babic will no longer depend 
on the coarse-graining length. The components of D, related 
to the characteristic size of the domain along each direction, 
are defined by the following equation: 



where Xi is the *-th coordinate of the averaging point. 

If the local homogeneity hypothesis holds, Di can be shown 
to scale with D m ^ 1 the size of the averaging domain in the 
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i direction (see Appendix B). The scaling will depend also 
on the shape of the weighting function. Moreover, if the 
assumption of local homogeneity holds, T 7 contains all the 
dependence on the coarse-graining width of the stress tensor 
while the kinetic stress tensor defined by means of the 
present fluctuation decomposition appears to be independent 
on averaging domain size. 

The translational kinetic energy balance is only slightly 
more complicated than the one obtained by Babic. Once 
more, most of the new terms are negligible under the local 
homogeneity assumption. 

The total translational kinetic energy is found to be the 
sum of three contributions, pEx + pE j, + pex, which are re¬ 
spectively the translational kinetic energy of the mean flow, 
the translational kinetic energy related to the velocity gradi¬ 
ents, and the truly fluctuational translational kinetic energy. 
These three variables are defined as 


pEj 1 = p^v • v, 

(13) 

ip(D • Vv) • (D • Vv), 

(14) 

POO 

/ WpUipVp ■ v p dt'. 

J —OO „ 

(15) 


The translational kinetic energy balance equation is derived 
as: 

d 

— (pE T + pE t + per) + V • {pE T + pE% + pet)v = 

V ■ (q k + q c ) 

+V • ((T fc + T c + T 7 ) • v) 

+pg • v - 7 fc - y 7 , ( 16 ) 

where q k and q c are Babic’s energy fluxes (see Appendix A). 
The last two terms correspond to the rate of conversion of 
kinetic energy to other forms of energy (dissipation due to 
friction, collisions, plasticity, but also to storage and resti¬ 
tution of elastic energy). In the rigid limit they represent 
dissipation through friction and inelastic collisions and are 
therefore positive. The first one contains the contribution of 
fluctuations: 

/ OO _ 

( v ~p ~ v ~q ) W pq dt 'i ( 17 ) 

■°° p q>p 

and the second one the rate of conversion of mechanical energy 
due to affine deformations: 

/ OO 

E E f P? ‘ (( x p ~ x u) ■ V v)Wp,dt', ( 18 ) 

■°° p q>p 

where w pq is the mean of w p and w q . Both terms are in¬ 
dependent of the size of the averaging domain. We can see 
here again the benefit of our approach: On one hand the do¬ 
main size-dependent terms are clearly identified in the equa¬ 
tion ( pE j, and T 7 ); on the other hand, the contribution of 
affine deformations is separated from the contribution of the 
fluctuations. 


In fact, by superficially interpreting Babic’s approach, one 
could conclude that all the energy is dissipated by means of 
fluctuations. This is incorrect, since affine deformations do 
also dissipate energy; the reason for this misunderstanding 
is that affine deformations are considered as fluctuations in 
Babic’s method. 

We can remove the mean flow kinetic energy terms by ma¬ 
nipulating the momentum equation in a classical way (for 
example, see Ref. Hi)- The result is an equation expressing 
the balance for the fluctuating energy related to fluctuations 
and to affine deformations: 

d 

ft (P E T + P £ t) + V • (pEj, + pe r)v = 

( T c + T 7 + T fc)t ; Vv + V • (q k + q c ) - 7 fe - y 7 , (19) 

where f stands for transposition. Let us do some more ma¬ 
nipulation. By looking at the term accounting for the rate of 
conversion of mechanical energy due to affine deformations, 
manipulation of Eq. (I18f) . making use of l p? = x 9 — x p , allows 
us to rewrite 7 7 in the following form: 

7 7 = T c * f : Vv, (20) 

where 

/ OO 

EE< Wpqdt', ( 21 ) 

■°° p q>p 

is a new contact stress tensor, differing from the T c by the 
weight used. Indeed, as mentioned above, w pq is the simple 
mean of w p and w q , whereas w pq is related to the fraction of 
the branch vector l pq lying within the averaging region. 

III. ONE DIRECTIONAL, STEADY, FULLY 
DEVELOPED FLOW 

In this section we will apply our method to a simplified flow 
situation in order to gather some information on the effect of 
the new terms on force balance and on energy cascade. 

Let us consider a one directional, steady, fully developed 
flow [defined by v = (y x , 0, 0), d x = 0, d t = 0[. In such a case, 
the continuity equation is identically zero. The momentum 
equation becomes: 

V • (T c + T 7 + T fc ) + pg = 0, (22) 

and the fluctuating energy equation 

(T° + T 7 + T fc )t : Vv + V • (q k + q c ) - q fc - q 7 = 0. (23) 

According to the definition of T 7 , it is easy to see that its 
only nonzero component is: 

T 7 , = - p(Dyd y v x + D z d z v x f , (24) 

and it is therefore straightforward to see that 

(T 7 ) t : Vv = 0, (25) 

that is, the contribution of affine deformations to the stress 
tensor does not perform any work. Moreover, it is obvious 
that T 7 directly affects only the x component of the force 
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balance. Given that the flow is fully developed, we can safely 
remove it out from the balances, obtaining: 

V ■ (T c + T fc ) + ,9g = 0, (26) 

( T c + T fe ) t : Vv + V • (q k + q c ) - 7 fe - 7 7 = 0. (27) 

According to Eq. (EU1) : 

(T c ) t : Vv — q 7 = (T c — T c *) : Vv. (28) 

It is clear that T c and T c * are not in principle the same, par¬ 
ticularly when the averaging domain is smaller than a particle 
diameter. However, if the averaging domain spans more than 
one particle diameter, we can assume that T c « T c *. Equa¬ 
tion m becomes therefore: 

T fc : Vv + V • (q k + q c ) - 7 fc = 0. (29) 

This means that the work provided by the contact part of 
the stress tensor is converted locally by affine deformations, 
while only the kinetic part of the stress tensor enters into a 
nonlocal energy balance. This result has important conse¬ 
quences, since it shows that separate paths may exist for the 
dissipation of the work of kinetic and contact stresses. How¬ 
ever, it does not cast doubt on the validity of hydrodynamic 
theories (Tij [73 since our results do not deal with constitutive 
relations but suggest that a revision of the balances used for 
such theories could lead to a simplified description in some 
cases. Concerning the size of the averaging domain, it may 
seem strange to talk about domains smaller than a particle 
diameter, because it may seem obvious that such domains 
are inappropriate: Few particles are indeed contained in a 
snapshot of such a domain. However, we are always deal¬ 
ing with space-time averages, so even if the averaging domain 
is small, in many cases (e.g., slowly evolving or steady-state 
processes) the domain indeed contains a lot of particles since 
many of them traveled through it during their history without 
affecting local homogeneity. In other cases, periodic bound¬ 
ary conditions may help too, allowing use of small domains 
along some directions. That is one of the reasons why in the 
next section we will employ the averaging method to some 
DEM simulation results pertaining to this class of simplified 
flows. 


IV. DEM SIMULATIONS 
A. Method & parameters 

We use our own 2D implementation [13 of the classical 
discrete element method where Newton’s equations of motion 
for a system of N “soft” disks are integrated. Such a 
technique, which is able to reproduce successfully the exper¬ 
imental results in many configurations (e.g., gravity-driven 
flows [Tl. 13. l2lLI 22 I]. sheared systems @,H3, granular materials 
close to jamming [24}, silos [25j, and rotating drums [26l [27| j. 
requires giving an explicit expression for the inter-particle 
forces. The discrete element method is classical and well 
known and can found in the aforementioned references. 
Therefore, we just present here the forces used in this work. 



FIG. 1: Snapshot of the simulated system. A 2D granular material 
is sheared between two bumpy walls. The lower wall is fixed, the 
upper one moves horizontally at a constant velocity Vt op = 100 
and is submitted to a constant vertical stress P = 1. 


For the normal force between two overlapping disks we 
use a standard linear spring-dashpot interaction model 
F n = k n S n — 7 nVraj where 8 n is the normal overlap, k n is 
the spring constant, 7n the damping coefficient, and v n the 
normal relative velocity. The damping models the dissipation 
characteristic of granular materials. Likewise, the tangential 
force is modeled as a linear elastic and linear dissipative 
force in the tangential direction F t = ktSt — 'JtVt, where k t 
is the tangential spring constant, St the tangential overlap, 
7 t the tangential damping, and Vt the tangential velocity 
at the contact point. The magnitude of St is truncated as 
necessary to satisfy the Coulomb law, |F t | < n\F n \, where /i 
is the grain-grain friction coefficient. The simulated system 
is two dimensional (Fig. [TJ) . The granular material is a 
dense assembly of N dissipative disks of average diameter 
d, and average mass m. A small polydispersity of ±20% is 
considered to prevent crystallization. 

The granular material is submitted to a plane shear, with¬ 
out gravity, leading to a uniform stress distribution. The ma¬ 
terial is sheared between two parallel rough walls, separated 
by a distance H. One of the walls is fixed, while the other 
moves at the prescribed velocity V. The flow and transverse 
directions are respectively called x and y. Periodic bound¬ 
ary conditions are applied along the flow direction, and the 
length of the simulation box, L x , is set to 60 grain sizes. 
The wall roughness is made of disks sharing the characteris¬ 
tics of the flowing grains (same polydispersity and mechanical 
properties, no rotation). Their centers are equally spaced by 
a distance equal to the largest disk diameter 1.2d. Along 
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(a) (b) 




FIG. 2: (Color online) (a) x velocity and (b) solid fraction profiles along the y direction ( Vt op = 100). The system appears to be in a 
simple shear condition (linear velocity profile, constant solid fraction) in the center of the cell. The size of the averaging window is 2d. 

(a) (b) 




FIG. 3: (Color online) (a) Fluctuating energy calculated by the present approach along the y direction ( Vt op = 100). The size of the 
averaging window is 2d. (b) Relative difference (in percentage) between the fluctuating energy profiles obtained in several ways [Babic’s 
(solid lines) or the present method (dashed lines), where the color represents the size of the averaging window] and that obtained with 
the present method and with an averaging window of Id. 


the y axis, the position y = 0 corresponds to the center of 
the glued grains on the fixed wall. The normal stress ap¬ 
plied on the moving wall, Poj is controlled. The vertical po¬ 
sition of the moving wall is thus not fixed and, using the 
method described in Ref. [28[, the height of the system H 
obeys H = (Pq — P w )L x /g p , where g p is a viscous damping 
parameter and P w is the normal stress exerted by the grains 
on the moving wall. 

The following values of the parameters are used: k n /Po = 
10 5 , kt = 2k n /7, 74 = 0, and y = 0.5 and g p = 100 \/k n m. 
The value of "f n is adjusted to obtain a normal restitution 
coefficient e n = 0.5 [5j. The equations of motion for the 
translational and rotational degrees of freedom are integrated 
with a velocity-Verlet scheme with a time step 10 \Jm/P q. 


The number of grains, N , is adjusted so the size of the system 
in the y direction, H , is roughly equal to 80d. The initial 
state of the system is a randomly diluted hexagonal lattice 
with disorder both on grain positions and velocities. The 
attainment of a steady state is verified by observing time- 
invariant total kinetic energy and distance between the two 
walls. 


B. Averaging procedure 

For the simulation set up described above, we implemented 
an averaging procedure following the method developed in the 
earlier sections of the paper. 
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Several snapshots of the state of the system (positions, ve¬ 
locities, contact forces) at given times were captured. For each 
time a space average was performed; given that the system 
was periodic, averages were computed at different y values. A 
Heaviside step function was chosen as the weighting function. 
The only parameter of the weighting function was therefore 
the averaging window size. Some authors suggested 
smoother weighting functions. Here we preferred a step func¬ 
tion for its simplicity and to avoid smoothing the profiles too 
much. Given that the system was simple and periodic, we felt 
that a nonsmooth function was the best choice, so as not to 
lose information about layering, localization, and so on. 

Given that the simulations displayed a stationary state, the 
time average was performed as a simple mean over the results 
obtained for each snapshot. Standard errors of the time aver¬ 
ages were estimated by the blocking method [2^ in order to 
take into account time correlation of the data. 

Note that, in order to calculate the full set of variables, two 
steps were necessary: a first step to calculate average values 
of solid fraction, velocity, and contact stresses and a second 
step to calculate terms containing velocity fluctuations. Two 
steps are needed since velocity fluctuations are defined with 
respect to the mean field. 


C. Results 

In the following we will apply the method to a simulation 
with a velocity of the top wall V top = 100d \JP/m ; as stated 
above, the flow is at constant pressure, so the vertical position 
of the top wall evolved freely to a stationary value, yt op « 
83.9 d. 

It is convenient to characterize such a shear using the so- 
called inertial number I that compares the typical time scale 
of microscopic rearrangements with the typical time scale of 
macroscopic deformations: I = 'yd,/yJP/ps- Note that / cor¬ 
responds to the square root of the Savage number [l4[ and is 
also called the Coulomb number [30 ]. Our system corresponds 
to a rapid flow with a global value of the inertial number cor¬ 
responding to I > 1; such a flow regime was chosen because 
kinetic and contact stresses have the same order of magni¬ 
tude, allowing us to test the results of the previous sections. 
A snapshot of the system at steady state is given in Fig. [T] 
All profiles shown in the following are calculated using an 
averaging window of 2d unless otherwise stated. 

Figure [3][a) shows the velocity profile at steady state. As 
highlighted in the figure, the profile is clearly linear except 
from two zones (~ 10d wide) near the bumpy walls, where 
the shear rate increases approaching the boundary. 

The solid fraction profile [Fig. [2])b)] displays a central zone 
with a constant value, and a decreasing behavior approaching 
the walls. The state of the system, characterized by a lower 
compaction near the walls, can be appreciated also from Fig. 
Ul In fact, all the variables suggest that the system is in 
a simple shear condition (uniform shear, stresses, no energy 
flux) in the center and deviates from this state only in two 
narrow bands (~ 10d wide) near the bumpy walls. Another 
feature of the system that can be appreciated in Fig. [2jb) is 
a slight asymmetry of the profiles near the walls. As it could 
be seen also in Fig. [T] in the shear zone near the upper wall 
the material is denser and less agitated than in the shear zone 


(a) 




(c) 



y/d 


FIG. 4: (Color online) Components of the stress tensor [contact 
part (a), kinetic part (b), and total tensor (c)] along the y direction. 
The upper wall velocity is Vtop = 100 and the size of the averaging 
window is 2d. 


near the bottom wall. This slight difference is related to the 
asymmetry of boundary conditions: while the bottom wall is 
fixed, the upper one is free to move vertically in order to im¬ 
pose a constant normal stress. This yields slightly asymmetric 
profiles for all the variables. 

Figure Ha) displays the fluctuating energy profile calcu¬ 
lated with the new procedure developed in this work for an 
averaging window of 2d. To compare the two approaches (the 
present one and Babic’s one) we report in Fig. HT) the quan¬ 
tity k — 1, where k is the ratio of the fluctuating energy profile 
calculated using one of the methods with different averaging 
sizes to a reference profile obtained with the present method 
and an averaging size equal to D = Id. It clearly shows that 
the Babic’s method leads to an estimate of the fluctuating 
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(a) (b) 




FIG. 5: (Color online) Kinetic (a) and contact (b) fluxes of fluctuating energy along the y direction (Vt op = 100). The size of the 
averaging window is 2d. 

(a) (b) 




FIG. 6: (Color online) Stress power and dissipative terms: (a) contact stress tensor and (b) kinetic stress tensor, along the y direction 
(Vtop = 100). The size of the averaging window is 2d. 


energy which increases with the domain size. On the con¬ 
trary, the present method leads to an estimate which does 
not depend on the size of the averaging window. However, if 
the latter size is too large it influences the averages close to 
the walls due to the inhomogeneity of the shear in such an 
averaging window. To illustrate this point, we have reported 
in Fig. [U)b) results obtained with D = 5 d. Note also that, 
using the present method, the fluctuating energy profiles for 
D = Id and D = 2d are almost equal for any position y. 

Figure |U)a) displays the profiles of the components of the 
contact part of the stress tensor (here mainly due to collisions, 
since the flow is quite rapid). Again, a plateau is found at 
the center; when approaching the walls all the components 
decrease in absolute value. 

At the same time [Fig. |U(b)], the components of the ki¬ 
netic part of the stress tensor increase when approaching the 
walls. This is a well-known direct consequence of the decreas¬ 
ing of solid fraction. In our simulation, near the walls, kinetic 
stresses are larger than collisional ones. 


The sum of the two parts gives the total stress tensor, shown 
in Fig. 0|c). As predicted by the stress balance, the (xy) and 
(yy) components of the total stress tensor are uniform along y. 
Concerning the third part of the stress tensor, T T , it is clear 
that according to Eq. (EH) its only nonzero component is the 
(. xx ) one. For the reference simulation it was verified that, 
apart for a small difference in a thin zone (~ 2d) close to the 
walls (where slight deviations from local homogeneity appear 
due to the shape of the wall), D y ss D m y /y/l2, confirming 
the scaling given above. 

Figure [5] collects contact and kinetic flux profiles. All the 
components of these fluxes differ from zero only when ap¬ 
proaching the wall; as regards the y component of the fluxes, 
fluxes are negative when the fluctuating energy decreases with 
y and positive otherwise. This agrees with the sign conven¬ 
tion used for q in the derivation, and supports a Fourier-like 
expression for the flux dependence on fluctuating energy gra¬ 
dient. On the other hand, the horizontal components of the 
fluxes are not negligible near the walls. This contrasts with a 






































(a) 


(b) 



FIG. 7: (Color online) Layering close to a bumpy wall: solid fraction computed using the present method with an averaging window of 
Id, on lines and on slices (Id wide) for (a) Vt op = 100, and (b) Vt op = 10. 


linear transport theory since there is no gradient of fluctuat¬ 
ing energy in the x direction. 

Figure [6] displays the estimates of the stress power related 
to contact and kinetic contributions and of the dissipative 
terms related to affine and non affine motions. Firstly, it can 
be seen that far from the walls, the algebraic relations hold 

T ct : Vv = T fct : Vv = 7 fc , (30) 

therefore justifying the separate dissipation paths for kinetic 
and contact stress powers. Then, near the walls, the energy 
fluxes seem to influence only the kinetic part, therefore sup¬ 
porting the result obtained in the previous section: only ki¬ 
netic terms enter in a nonlocal energy balance, while contact 
stress power is dissipated locally. 

V. DISCUSSION 

In the previous sections we developed an averaging proce¬ 
dure for granular materials which was based on a weighted 
average plus an original fluctuation decomposition. Through 
numerical simulation results we proved the strength of the 
approach which gives consistent predictions for kinetic terms, 
and—more importantly—gives a new picture of the energy 
cascade in granular flows, at least for simple flows. In this 
section we discuss the validity of such an approach. It is clear 
that the assumption of local homogeneity holds only when the 
solid fraction and the shear rate profiles do not display large 
gradients with respect to the averaging domain size. For the 
simulations presented above, two zones were identified near 
the walls where such gradients exist. The independence of the 
fluctuating energy profile on the size of the averaging window 
and the respect of the stress balance are two proofs that the 
assumption is still acceptable. 


Figure [7] allows us to discuss another important point: the 
physical meaning of the variables obtained by the method, 
with a particular regard to solid fraction. In the figure, the 
solid fraction profile was computed by two other methods: 
estimating the fraction of (1) a horizontal line and (2) of a 
horizontal slice (Id wide) intersecting the particles. These 
two methods have a clear geometrical meaning, the second 
being a moving average of the first. 

For the rapid flow simulation which was taken as the refer¬ 
ence [Fig. [TJa)], where no layering occurs, the three methods 
agree very well. FigurtJTKb) displays the results of the three 
methods for a denser flow (V top = 10, y top ~ 50, => I ~ 0.2). 
Please note that the quite ordered bumpy wall was explic¬ 
itly prepared to display layering. In this case the estimates 
on lines strongly fluctuate near the wall, and the average on 
slices is a smoother profile which still fluctuates, while the 
estimate produced by the present averaging method seems 
to be intermediate between the two. Far from the wall re¬ 
sults from the methods coincide. Therefore, in the presence 
of layering, the solid fraction calculated by use of the present 
method may not have the direct geometrical meaning that its 
name implies. This problem may be reduced by choosing a 
smoother weighting function (in this paper we chose on pur¬ 
pose a very abrupt weighting function); we also suspect that 
in three dimensions this rather pathological issue may be less 
important. 

Another issue caused by layering is that the strong spatial 
fluctuation of variables may imply that some of the assump¬ 
tions of the local homogeneity hypothesis do not hold: The 
distribution of particles around the averaging point is not in 
principle homogeneous. Fluctuations are evidently associated 
with strong local gradients of solid fraction and shear rate. 
Therefore, when strong layering occurs, care has to be taken 
in applying the present method due to the approximations 
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therein. 

The inhomogeneous distribution of particles inside the av¬ 
eraging domain yields also a deviation of the virtual center 
of mass X from the center of the domain. This is not, in 
principle, a problem, because the averaging procedure is re¬ 
lated to a volume and not to a point; average variables may 
be harmlessly referred to the virtual center of mass instead 
of the simple geometric center of the averaging domain. In 
other words the relation X = x which seemed to be a conse¬ 
quence of the local homogeneity assumption can be taken as 
a definition of x, the point where averages are supposed to be 
related. 

In addition, results from simulations at V top = 1 and V top = 
10 (which show quite strong layering), not shown here for 
the sake of brevity, show that stress balance is respected and 
confirm locality of the dissipation of the contact stress power. 
However, if layering occurs, local homogeneity may be difficult 
to assume, and therefore the simplified balances obtained in 
this work have always to be critically analyzed when dealing 
with regions very close to the walls. 


VI. CONCLUSIONS 


In this work, we propose a new derivation of continuum bal¬ 
ance equations for granular materials which is inspired by the 
work of Babic |7|. We introduce a new decomposition of par¬ 
ticle scale velocities into mean and fluctuations, taking into 
account velocity gradients. A very reasonable local homogene¬ 
ity hypothesis allows us to simplify the balances, ensuring at 
the same time representativity and physical meaning of the 
continuum estimates obtained. 

Two important results are obtained. First, the method 
solves the problem of dependence of some variables (e.g., fluc¬ 
tuating energy, kinetic stresses) on the averaging domain size 
which is intrinsic to previous methods Q by clearly identifying 
scale-dependent terms. Second, the development suggests a 
new physical picture for the energy cascade in granular flows: 
kinetic terms enter in a nonlocal energy balance, while con¬ 
tact terms (containing both collisions and long-lasting con¬ 
tacts) seem to be dissipated locally. Both results are verified 
against discrete element simulations of granular flows, with 
respect to which the method is shown at work. 

This point is particularly important with regards to the 
hope of obtaining a statistical mechanics for granular flows 
whose predictions can be faithfully tested in physical experi¬ 
ments and numerical simulations. Note also that the question 
of the size of the averaging window close to a boundary [TlJ is 
problematic since it is bounded by the distance between the 
center of a grain and that boundary. 

A last note has to be added on the domain of application of 
the present method. Discrete element methods are sometimes 
divided into smooth and nonsmooth ones a The present 
averaging method (as Babic’s one) applies rigorously only to 
smooth methods, since the idea of smooth (differentiable) par¬ 
ticle fields was implicit in the expression of the equation of 
motion. The adaptation of such a method to non smooth dy¬ 
namics is a very important point which will be the subject of 
future works. 
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Appendix A: Derivation of balance equations 

The derivation follows closely that of Babic, except for the 
choice of the fluctuation decomposition. The evolution equa¬ 
tion for a generic particle property i/j p is: 


d4\ 


a rp _ \ ' p 
dt ~ ^ p ' 


m p g p . 


We multiply by w p , sum on p, and integrate on time: 

r sr ^ 

/ 2^ w P m p~dir dt = 


(Al) 


(A2) 


EE w p P pq dt' + / ^ w p m p g p dt'. 

J-oo p q J -oo p 

The left hand side term can be developed into 0: 

d f °° . r°° 

AT / E w P m P i’ P dt' + V • / Y w p m p ^pdt', (A3) 

J —oo „ J — oo „ 


As regards the right hand side of Eq. (|A2f) . Babic also 
showed that the interaction term could be developed as: 


E E w p^pqdt — 


„ i r 
+V '2/_„ 


~ p q 

/•OO 

/ ^ ' 'y ' (Ppq A Pqp) w pqdt 

— p q>p 

E E W (Ppq ~ Pqp) w pqdt ■ (A4) 

p q>p 


where w pq and w pq are two weighting functions respectively 
defined as: 


■Wpq = 


/ wj(x p + sip, — x, t' — t)ds, 

Jo 


and 


g _ W p + W q 


(AS) 


(A6) 


Now let us see the effect of the fluctuation decomposition 
on Eq. (IA3I) . Decomposing the velocity field with the new 
decomposition leads to: 

Vp(xp) = v(x) + (Xp - x) • Vv(x) + Vp(xp). (A7) 
The divergence term on Eq. (IA3I) can be written as: 

V • p\rJ[) + V • pi/iv 

TV ' (^J E w p m p' l l ) p{*-p ~ ■ Vv^ (A8) 
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Decomposing also ip p yields for the second term in Eq. (1A8[) : 

pipv = I / w p m p Vp(xp — x) ) • V^(x)dt' 

V~°° p ) 

/ OO 

X Wpm p Vpijjpdt. (A9) 

-OO ^ 


Given that 


X*, w p m p (x p — x)dt / = p(X — x), 


(A10) 


where X is the average center of mass of the averaging domain, 
the integral in the third term of Eq. (IA8I) becomes 


/ OO 

X w p m p ’i/jp(xp — x)dt' = -i/5p(X — x) 

-OO p 

/ OO 

y^rc p mp((xp — x) • V'0)(xp — x)df' 

-oo p 

/ OO 

X, w p m p ip p (xp - x.)dt r . (All) 

-OO ^ 


Following the homogeneity assumption, we neglect correla¬ 
tions of fluctuations and of particle positions, and we consider 
that the average center of mass coincides with the averaging 
point. This means that terms under the divergence operator 
in Eq. (IA8I) simplify to 



WpTOp((x p 


- x) • V'0)(Xp 


pvi/j + pipV + 

— x)dt' ) • Vv. (A12) 


Developing the k component of the third term in Eq. (IA12I) 
[using the definition Ax p i = ( x p i — xp)] yields: 


/ OO 

wj p mp((x p - x) • V^)((xp - x) • S7v k )dt' = 

-OO ^ 


/ OO 

X - 

-oo ^ 


WpmpAxpiAxpjdt' .(A13) 


* j 


Assumption of local homogeneity also implies decorrelation 
of the components of particle position with respect to the 
averaging point. This means that: 

/ OO 

^WpmpAxpiAxpjdt' = pD^ if i — j, 0 otherwise, 

-oo p 

(A14) 

where Di is related to the scale of the averaging volume in 
the i direction. If D is the vector with components Di, the 
third term of Eq. (IA12I) is 

p(D • V^)(D • Vufc). (A15) 


This term renders explicit the scale dependence which was 
otherwise implicit in Babic’s method. Equation <1A3ll is finally 
reduced to: 

d - 

—pip + V • pvip + V • pipv + V • [p(D ■ VV’)(D ■ Vv)] . (A16) 


As regards to the effect of the fluctuation decomposition on 
the right hand side of Eq. (1A2[) . the result depends on the 
particular property ipp considered. Thus in the following we 
will analyze each property separately. 

If mass conservation is considered (ip p = 1, P pq = 0, 
g p = 0), the right-hand-side term of Eq. (1A2I) is zero and it is 
straightforward to see that the classical continuity equation 
is obtained [Eq. ([5]. 

If force balances are considered (ip p = v p , Ppq = f pq , g P = 
g), Eq. (IA3I) becomes 

d _ 

—pv -f V ■ pvv + V ■ pvv + V • [p(D ■ Vv)(D • Vv)l. (A17) 

at 

The last two terms can be considered as related to two com¬ 
ponents of an effective stress tensor as defined by Eqs. m 
and m- On the other hand, the right-hand side of Eq. (IA2I) 
is not affected by the nature of the fluctuation decomposition. 
Substitution of the P pq = f pq and g p = g yields to the stan¬ 
dard momentum equation, Eq. ©, with the third component 
of the effective stress tensor given by Eq. ©• 

Let us derive now the translational kinetic energy equation. 
In this case the reference particle property is 

= \v p ■ v p . (A18) 

The average translational kinetic energy is 

pip = 7) Yw p m p -v p • v p dt', (A19) 

Z p 

and the fluctuation decomposition allows us to split this quan¬ 
tity into six terms, three of which disappear due to the local 
homogeneity assumption, yielding: 

ipv-v 

+ 9 / X w p m p(( x p - x ) • V ^) • (( X P - x ) • Vv)dt' 

ZJ -oo p 

1 f°° 

+ 9 / X ■ Vprff' = pE T + pEPp + pe T (A20) 

1 p 

where, with respect to Babic’s development, per is likely 
to be scale independent, and a new term appears, which is 
the translational kinetic energy related to affine deformation 
within the average volume pEj. 

The divergence term in Eq. (IA3I) corresponds, for the ki¬ 
netic energy, to: 

/ oo i 

X “’p m p^( v p ' Vp)v p dt'. (A21) 

-oo p 1 

When developing this term, the fluctuation decomposition 
gives 27 terms, most of which can be deleted as in the pre¬ 
vious development through (1) the postulate X = x, (2) the 
assumption of decorrelation among components of x p — x, 
and (3) the assumption of decorrelation between components 
of Xp — x and velocity fluctuations. It can be shown that the 
irreducible terms are 

(pE T + pElp + /9£t)v - q k - T fc • v - T 7 • v, 


(A22) 
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where 


which correspond respectively to: 


q 


k 


1 f°° 

2 E w p m p (y p ■ v p )v p df'. 


(A23) 


Let us treat the terms coming from the right hand side 
term of Eq. (IA2I) . First, we must identify the particle scale 
interaction terms. Taking the product of v p and Newton’s 
equation, 


UVp , \ -\ . 

to p v p • — = v p ■ CL f p<? + TO p g ) 

9 

=> Vp = v p • (X! f pq + m p g )- ( A24 ) 

9 

Therefore, the particle interaction term is 

El fp9 • Vpi (A25) 

9 

and the particle source term is 

m p g-Vp. (A26) 

Concerning the source term, local homogeneity implies: 

/ OO 

E w p m p g ■ v p dt' = pg ■ v + (p(X - x) • Vv) g = pg-v. 

(A27) 

Developing the interaction terms yields: 


- 7 fe - 7 7 + V-q c + V-(T c -v), (A30) 

where the hrst two terms are the conversion rates of transla¬ 
tional kinetic energy into other forms of energy respectively 
due to velocity fluctuations (— 7 fc ) and to affine deformations 
(— 7 7 ); q c can be seen as a translational kinetic energy flux 
related to contact forces. 

The translational kinetic energy balance equation is there¬ 
fore given by Eq. m, which can be further simplified to 
Eq. (1191) when subtracting the mean flow energy balance. 

Appendix B: Derivation of a scaling for D 

In the following, a simplified scaling for the vector D which 
enters in the definition of T 7 is derived using the local homo¬ 
geneity hypothesis. Let us assume equal-sized spheres, and a 
step function on space and time as a weighting function. Let 
us consider variables computed at a point (x, t ) in space and 
time. The average density is therefore: 


/ OO 

'22’w p dt', (Bl) 

-OO p 

where w p = 1/(VT) for particles residing in the averaging 
volume V during the period (f — T/2, f + T/2) and zero other¬ 
wise. If we call N p the time average of the number of particles 
lying in the averaging volume, it is clear that p = m p N p /V. 
Developing the definition of the vector D, we find: 


EE( f P « • Vp + f qp • V q ) W pq dt' 

V 9 >P 

+v • j I y ] y 1 (fp 9 ' v p — ' v q) ^pq w P qdt ■ (A28) 

“ d-oo p q>p 

Using the fluctuation decomposition, the relation f pq = — i qp , 
and the assumption of local homogeneity this becomes: 


pDf = m p Ywpjxpj - Xj) 1 2 3 4 dt', (B2) 

J-oo p 

where x p t is the i component of the particle position. The 
local homogeneity hypothesis says that X = x and that x P i 
is uniformly distributed around 27 . This yields for large N p T 
to: 




E E ^9 ' ( v p 

p 9 >p 


)w° q dt' 


+ 



EEw (( x p 

P 9 >P 


Xq) • S7v)w pq dt' 




EE fp 9 * (Vp T Vq)l p qW p qdt 

P 9 >P 


+V- 



E E ^pq^pq w P qdt 

p q>p 



Therefore if the local homogeneity holds, Di is likely to scale 
on the size of the averaging domain in the i direction. The 
proportionality will then depend on the choice of the weight- 
(A29) ing function, on polydispersity, on shape, and so on. 
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